##########################################
# Aligned adjacent municipality analysis
##########################################

rm(list=ls())

#sink("~/Dropbox/Crime Chile/11_replication/06_spillover_crime.txt")

library(Hmisc)
library(ggplot2)
library(stargazer)
library(foreign)
library(rdrobust)
library(rdd)
library(readstata13)
library(gridExtra)

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<40){
    warning("Fewer than 40 clusters, variances may be unreliable")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

################
# Prepare data 
################

# reada data
d = read.dta13("~/Dropbox/Crime Chile/11_replication/local_crime_data_chile_2022january.dta")   
names(d)

# control
d$t_spill2 = 1 - d$t_spill

# spill
table(d$t_spill)

table(d$total_adjacent_county)

##############################################################################################
# Aligned adjacent municipalitys 1: Interaction with binary indicator: Property street crime
##############################################################################################

# main regression
d2 = d[!is.na(d$propstreetcrime_s),]
d2 = d2[d2$total_adjacent_county!=0,]
h1 = rdbwselect(d2$propstreetcrime_s,d2$margin,cluster = d2$cluster)$bws[1]
X = d2$margin
center = 0
bw = h1
d2$kweights = kernelwts(X, center, bw, kernel = "triangular")
sum(d2$kweights!=0)
h1

# interaction 
reg1 = lm(d2$propstreetcrime_s ~ d2$ruling + d2$t_spill + d2$t_spill*d2$ruling + d2$margin, weights = d2$kweights)
reg1_c <- round(coeftest(reg1, vcov = vcovCluster(reg1, cluster = d2$cluster)),4)
reg1_c 

stargazer(reg1_c, 
          summary=FALSE, 
          rownames=FALSE, 
          omit = c("Constant", "margin"),
          covariate.labels = 
          c("Political alignment (beta 1)", "Aligned adjacent municipality (beta 2)", "Interaction (beta 3)"),
          add.lines = list(c("MSE bandwidth","0.093"),c("N","665")),
          out="~/Dropbox/Crime Chile/11_replication/tableA13.html")

sink()
